/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */



/*! \file delta_particle.h
    \brief Multidimensional generalization of particle<OS_double> Laplacian operators.

	Multidimensional generalization of particle<OS_double> Laplacian operators
	\f[
			\delta_{S} u_i= \sum_{j \in \nu_{i,j}} \frac{S_i-S_j}{|S_i-S_j|^2}
	\f]
*/
#if !defined(_ParticlesDelta_)
#define _ParticlesDelta_
#include <string>
#include <src/miscellaneous/OS_basic.h>
#include <src/math/Operators/ParticlesOperators.h>
#include <src/link/Qhull_link.h>
#include <src/math/BoundaryCondition/ClassicalGraphDefinition.h>

#pragma warning( disable : 4355)
/*! 
    \brief General case for particle<OS_double> Laplacian operators.

	Multidimensional generalization of particle<OS_double> Laplacian operators
	\f[
			\delta_{S} u_i= \sum_{j \in \nu_{i,j}} \frac{S_i-S_j}{|S_i-S_j|^2}
	\f]
*/
template <class data, std::size_t = data::Dim>
class Delta_p {
		OS_STATIC_CHECK(false); // should never be instanciated
};

template <class data>
class Delta_p<data,2> : 
	public QHULL::qh_edges_set_type<data,data>
{
public :
		typedef QHULL::qh_edges_set_type<data,data>						Base;
		typedef  Particle_Operator										Family_Operator;
		typedef  typename Base::Base_Cast								Base_Cast;
		typedef Delta_p													type;
		typedef typename Base::graph_type								graph_type;

		Delta_p( ) {};
		Delta_p(const Delta_p & data) : Base(data) {};
		Delta_p(const Base_Cast &first) : Base(first)
		{};
		virtual ~Delta_p(){};
		static std::string get_name(){return " Delta_{S}u ";};
		Delta_p & operator = (const Delta_p & from) {
			Base::operator =(from);
			return *this;
		};

		virtual smart_ptr_type Create_Instance() const {
			return Create_Instance_helper_arg<type,weight_set_type>(&weight_set_);
		};
		virtual void update(){
			Base::update();
			construct<weight_set_type>();
		};
		virtual const graph_type& get_graph() const {
			return graph_;
		};

		template <class T>
		void construct()
		{
			OS_STATIC_CHECK(false); // you must implement the Delaunay graph for your own weight set type
		};
		template <>
		void construct<weight_set_default_def<base_value_type>>() {
			typedef weight_set_type::iterator1 iterator1;
			typedef weight_set_type::iterator2 iterator2;

			graph_.clear();



			OS_int numpoints = get_state_variable()->size(2);
			base_value_type dim = get_state_variable()->size(1);
			cardinal_.resize(numpoints);

			std::fill(cardinal_.begin(),cardinal_.end(),0);

			iterator1 itl;
			iterator2 itl_;
			iterator1 end(weight_set_.end1());
			for (itl = weight_set_.begin1();itl != end;++itl) {
					size_type size = itl.index1();
					iterator2 end2(itl.end());
					for (itl_ = itl.begin();itl_ != end2;++itl_) {
							++cardinal_[size];
					};
					
			};

			for (itl = weight_set_.begin1();itl != end;++itl) {
					indice_type i = itl.index1();
					base_value_type size_i = cardinal_[i];
					iterator2 end2(itl.end());
					for (itl_ = itl.begin();itl_ != end2;++itl_) {
							indice_type j = itl_.index2();
							base_value_type size_j = cardinal_[j];
							(*itl_) = dim *( .5 / size_i + .5 / size_j);
							graph_.insert(std::make_pair(i,j));
					};
					
			};
		};

		std::vector<size_type> cardinal_;
};
/*! 
    \brief One dimensional specialisation of particle<OS_double> Laplacian operators.

	Compute 
	\f[
			\delta_{S} u_i= \frac{1}{S_i-S_{i+1}}+\frac{1}{S_i-S_{i-1}}
	\f]
*/

template <class A>
class Delta_p<A,1>
	: public QHULL::qh_edges_set_type<A,A>
{
public :
		enum:std::size_t												{Order = 2};
		typedef OTSCentered												Centered_Type;
		typedef Dirichlet												Boundary_type;
		
		typedef QHULL::qh_edges_set_type<A,A>							Base;
		typedef  typename Base::Base_Cast								Base_Cast;
		typedef Delta_p											type;

		Delta_p( ) : size_(0) {};
		Delta_p(const Delta_p & data) : size_(data.size_),  Base(data) {};
		virtual ~Delta_p(){};
		static std::string get_name(){return " Delta_{S}u ";};
		Delta_p(const Base_Cast &first) : size_(0), Base(first)
		{};
		Delta_p & operator = (const Delta_p & from) {
			Base::operator =(from);
			return *this;
		};
		virtual result_type operator()(typename argument_type from) const 
		{
			return ( !( from-(from>>1) ) + !( from-(from<<1) ) );
		};
		virtual void update() {
			Base::update();
			state_variable_type s_v = get_state_variable();
			if (s_v->size() != size_) {
				size_ = s_v->size();
				if (size_ > 0)
					construct(weight_set_);
					build_graph();
			};
		};
		template <class T>
		void construct(T&) {
			OS_STATIC_CHECK(false);//("you have to implement the rebuild function for your own weight set type");
		};
		template <>
		void construct(weight_set_default_def<base_value_type> & weight_set) {
			weight_set.resize(size_,size_,false);
			for (OS_size i = 0 ; i < size_ ; ++i) {
					indice_type k = (indice_type) (i-1);
					if ( k >=0 )
						weight_set(i,k) = 1.;
					k = i+1;
					if ( k < size_ )
						weight_set(i,k) = 1.;
			};
		};

		virtual smart_ptr_type Create_Instance() const {
			return Create_Instance_helper<type>();
		};
protected :
		size_type size_;
};

/*! \brief Common Bind Time operator for Laplacian particle<OS_double> operator Delta_p

	This is the time Bind operator of Laplacian difference operator Delta_traits, that is Delta_t(t,u)=Delta(u)
*/
template <class data, typename Time = OS_double >
struct Delta_p_t : public BindTimeOperator<Delta_p<data, data::Dim> ,Time > {

	typedef Delta_p_t														time_operator_type;
	typedef BindTimeOperator<time_operator_type ,Time >						Base;
	Delta_p_t(){};
};


#endif